#ifndef __Tempest_MATRIX__
#define __Tempest_MATRIX__

#include <iostream>
#include <new>
#include <vector>
#include <exception>
#include <cmath>
#include <limits>

/// Set the machine epsilon.
const double eps = std::numeric_limits<double>::epsilon();

template <typename T>

class Matrix
{
//构造矩阵类，也作为行向量类和列向量类这两个派生类的基类
private:
    int n_row,n_col;/*分别定义行数和列数*/
    T* matrix_elements=nullptr;/*定义一个一维指针，指向数组首地址，该数组用于存放矩阵各元素的值*/
    
public:
    Matrix(){};
    Matrix(int row,int col):n_row(row),n_col(col)/*构造函数，产生一个行数n_row，列数n_col的矩阵*/
    {
        if(n_row<1||n_col<1)
        Matrix<T>::warning(1);
        else
        {
            matrix_elements=new T [n_row*n_col];
        }
    };
    Matrix(const Matrix<T>& M):n_row(M.n_row),n_col(M.n_col)/*构造函数，复制功能*/
    {
        matrix_elements=new T [n_row*n_col];
        if(n_row==0||n_col==0)
            warning(2);
        else
            for(int i=0;i<n_row*n_col;i++)
                matrix_elements[i]=M.matrix_elements[i];
    };
    void Matrix_transport();/*矩阵转置功能*/
    bool is_square();/*判断是否为方阵功能*/
    void printMatrix();/*打印矩阵功能*/
    friend std::ostream& operator<<(std::ostream& out,Matrix<T>& M);/*重载流输出函数*/
    friend std::istream& operator>>(std::istream& in,Matrix<T>& M);/*重载流输入函数*/
    T count_rows()const{return n_row;};/*给出矩阵的行数*/
    T count_cols()const{return n_col;};/*给出矩阵的列数*/
    void reset(int row,int col,const T& _val);/*修改矩阵某个位置的值*/
    void makeEmpty();/*使矩阵各个元素归零，矩阵行列数不变*/
    void releaseMemory();/*释放内存*/
    Matrix<T>& operator+(const Matrix<T>& M);/*重载加法运算符*/
    Matrix<T>& operator-(const Matrix<T>& M);/*重载减法运算符*/
    Matrix<T> operator-();/*求负矩阵*/
    Matrix<T> Matrix_inv(const Matrix<T>& M);/*求矩阵的逆*/
    bool operator==(const Matrix<T>& M);/*重载==运算符*/
    bool operator!=(const Matrix<T>& M);/*重载！=运算符*/
    Matrix<T> operator*(const Matrix<T>& M);/*重载乘法运算符*/
    T &operator()(int _i,int _j);/*重载函数调用运算符*/
    void warning(int _error);/*警告功能，用于提示是否越界或者矩阵运算中是否出现size不匹配的情况*/
};

template <typename T1,typename T2>
class RowVector : public Matrix<T1>
{
//继承派生类，行向量
private:
    int n_row=1,n_col;
    T1* matrix_elements=nullptr;
public:
    RowVector<T1,T2>& operator*(const RowVector<T1,T2>& M);
};

template <typename T1,typename T3>
class ColVector : public Matrix<T1>
{
//继承派生类，列向量
private:
    int n_row,n_col=1;
    T1* matrix_elements=nullptr;
public:
    ColVector<T1,T3>& operator*(const ColVector<T1,T3>& M);
};

template <typename T>
void Matrix<T>::warning(int _error)
{
    switch(_error)
    {
    case 1:
    std::cerr<<"Out of range!"<<std::endl;
        break;
    case 2:
    std::cerr<<"Unequal size!"<<std::endl;
    default :
    std::cerr<<"Unexpected error."<<std::endl;
    break;
    }
}

template <typename T>
void Matrix<T>::Matrix_transport()
{
    int i,j;
    T* M;
    M=new T*[n_row*n_col];
    for(i=0;i<n_row*n_col;i++)
        M[i]=matrix_elements[i];
    delete [] matrix_elements;
    int change=n_row;
    n_row=n_col;
    n_col=change;
    matrix_elements=new T*[n_row*n_col];
    for(i=0;i<n_row;i++)
    for(j=0;j<n_col;j++)
        matrix_elements[i*n_col+j]=M[j*n_row+i];
    delete [] M;
}

template <typename T>
bool Matrix<T>::is_square()
{
    if(n_row==n_col)
        return true;
    else
        return false;
}

template <typename T>
void Matrix<T>::printMatrix()
{
    if(n_row==0||n_col==0)
    return;
    for(int i=0;i<n_row;i++)
    {
        for(int j=0;j<n_col-1;j++)
        std::cout<<matrix_elements[i*n_col+j]<<" ";
    std::cout<<matrix_elements[i*n_col+n_col-1]<<std::endl;
    }
}

template <typename T>
T &Matrix<T>::operator()(int _i,int _j)
{
    if(_i<0||_j<0||_i>=n_row||_j>=n_col)
        warning(1);
    return matrix_elements[_i*n_col+_j];
}

template <typename T>
std::ostream& operator<<(std::ostream& out,Matrix<T>& M)
{
    for(int i=0;i<M.count_rows;i++)
    {
    for(int j=0;j<M.count_cols;j++)
        out<<M(i,j)<<" ";
    out<<std::endl;
    }
    out<<std::endl;
    return out;
}

template <typename T>
std::istream& operator>>(std::istream& in,Matrix<T>& M)
{
    for(int i=0;i<M.count_rows;i++)
    {
    for(int j=0;j<M.count_cols;j++)
        in>>M(i,j);
    }
    return in;
}

template <typename T>
void Matrix<T>::reset(int row,int col,const T& _val)
{
    if(row<0||col<0||row>=n_row||col>=n_col)
    {
    warning(1);
        return;
    }
    matrix_elements[row*n_col+col]=_val;
}

template <typename T>
void Matrix<T>::releaseMemory()
{
    if(matrix_elements!=nullptr)
        delete [] matrix_elements;
    matrix_elements=nullptr;
    n_row=n_col=0;
}

template <typename T>
Matrix<T>& Matrix<T>::operator+(const Matrix<T>& M)
{
    if(n_row!=M.count_rows||n_col!=M.count_cols)
    warning(2);
    Matrix<T> temp(*this);
    if(n_row==M.count_rows&&n_col==M.count_cols)
    {
    for(int i=0;i<n_row*n_col;i++)
        temp.matrix_elements[i]=matrix_elements[i]+M.matrix_elements[i];
    return temp;
    }
    else
    warning(2);
}

template <typename T>
Matrix<T>& Matrix<T>::operator-(const Matrix<T>& M)
{
    if(n_row!=M.count_rows||n_col!=M.count_cols)
    warning(2);
    Matrix<T> temp(*this);
    if(n_row==M.count_rows&&n_col==M.count_cols)
    {
    for(int i=0;i<n_row*n_col;i++)
        temp.matrix_elements[i]=matrix_elements[i]-M.matrix_elements[i];
    return temp;
    }
    else
    warning(2);
}

template <typename T>
Matrix<T> Matrix<T>::Matrix_inv(const Matrix<T>& M)
{
    Matrix<T> temp(M);
    Matrix<T> eye(M.n_row,M.n_col);
    double ch;
    for(int i=0;i<M.n_row;i++){
        eye.reset(i,i,1.0);
    }
    for (int i=0;i<=M.n_row-1;i++)
        {
            int k;
            for (k=i;k<=M.n_row-1;k++)
            {
                if (fabs(temp(k,i))>1e-10)
                    break;
            }
            if (k <= M.n_row-1)
            {
                if (k != i)
                {
                    for (int j = 0; j <= M.n_row-1; j++)
                    {
                        ch = temp(i,j);
                        temp.reset(i, j, temp(k,j));
                        temp.reset(k, j, ch);
                        ch = eye(i,j);
                        eye.reset(i, j, eye(k,j));
                        eye.reset(k, j, ch);
                    }
                }
                double b = temp(i,i);
                for (int j = 0; j <= M.n_row-1; j++)
                {
                    temp.reset(i, j, temp(i,j)/b);
                    eye.reset(i, j, eye(i,j)/b);
                }
                for (int j = i + 1; j <= M.n_row-1; j++)
                {
                    b = -temp(j,i);
                    for (k = 0; k <= M.n_row-1; k++)
                    {
                        temp.reset(j, k, temp(j,k)+b*temp(i,k));
                        eye.reset(j, k, eye(j,k)+b*eye(i,k));
                    }
                }
            }
        }
        for (int i = M.n_row-1; i > 0; i--)
        {
            for (int j = i - 1; j >= 0; j--)
            {
                double b = -temp(j,i);
                temp.reset(j, i, 0);
                for (int k = 0; k <= M.n_row-1; k++)
                {
                    eye.reset(j, k, eye(j,k)+b*eye(i,k));
                }
            }
        }
    return eye;
}

template <typename T>
Matrix<T> Matrix<T>::operator-()
{
    Matrix<T> temp(*this);
    for(int i=0;i<n_row*n_col;i++)
        temp.matrix_elements[i]=matrix_elements[i]*(-1);
    return temp;
}

template <typename T>
bool Matrix<T>::operator==(const Matrix<T>& M)
{
    bool equal=true;
    if(n_row!=M.count_rows||n_col!=M.count_cols)
    warning(2);
    if(n_row==M.count_rows&&n_col==M.count_cols)
    {
    for(int i=0;i<n_row*n_col;i++)
        if(matrix_elements[i]!=M.matrix_elements[i])
            equal=false;
    return equal;
    }
    else
    warning(2);
    return false;
}

template <typename T>
bool Matrix<T>::operator!=(const Matrix<T>& M)
{
    bool equal=false;
    if(n_row!=M.count_rows||n_col!=M.count_cols)
    warning(2);
    if(n_row==M.count_rows&&n_col==M.count_cols)
    {
    for(int i=0;i<n_row*n_col;i++)
        if(matrix_elements[i]!=M.matrix_elements[i])
            equal=true;
    return equal;
    }
    else
    warning(2);
    return true;
}

template <typename T>
Matrix<T> Matrix<T>::operator*(const Matrix<T>& M)
{
    if(n_col!=M.n_row)
        warning(2);
    Matrix<T> temp(n_row,M.n_col);
    for(int i=0;i<n_row;i++)
        for(int j=0;j<M.n_col;j++)
        {
            temp(i,j)=0;
            for(int k=0;k<M.n_row;k++)
                temp(i,j)+=matrix_elements[i*n_col+k]*M.matrix_elements[k*M.n_col+j];
        }
    return temp;
}

template <typename T1,typename T2>
RowVector<T1,T2>& RowVector<T1,T2>::operator*(const RowVector<T1,T2>& M)
{
    if(n_col!=M.n_col)
        std::cerr<<"Unequal size!"<<std::endl;
    RowVector<T1,T2> neiji=0;
    for(int i=0;i<n_col;i++)
        neiji+=matrix_elements[i]*M.matrix_elements[i];
    return neiji;
}

template <typename T1,typename T3>
ColVector<T1,T3>& ColVector<T1,T3>::operator*(const ColVector<T1,T3>& M)
{
    if(n_row!=M.n_row)
        std::cerr<<"Unequal size!"<<std::endl;
    RowVector<T1,T3> neiji=0;
    for(int i=0;i<n_row;i++)
        neiji+=matrix_elements[i]*M.matrix_elements[i];
    return neiji;
}

#else
/// DO NOTHING.
#endif
